Minimizing total weighted latency in home healthcare routing and scheduling with patient prioritization

We study a home healthcare routing and scheduling problem, where multiple healthcare service provider teams should visit a given set of patients at their homes. The problem involves assigning each patient to a team and generating the routes of the teams such that each patient is visited once. When patients are prioritized according to the severity of their condition or their service urgency, the problem minimizes the total weighted waiting time of the patients, where the weights represent the triage levels. In this form, the problem generalizes the multiple traveling repairman problem. To obtain optimal solutions for small to moderate-size instances, we propose a level-based integer programming (IP) model on a transformed input network. To solve larger instances, we develop a metaheuristic algorithm that relies on a customized saving procedure and a general variable neighborhood search algorithm. We evaluate the IP model and the metaheuristic on various small-, medium- and large-sized instances coming from the vehicle routing literature. While the IP model finds the optimal solutions to all the small- and medium-sized instances within three hours of run time, the metaheuristic algorithm achieves the optimal solutions to all instances within merely a few seconds. We also provide a case study involving Covid-19 patients in a district of Istanbul and derive insights for the planners by means of several analyses.


Introduction
Home healthcare services have been growing globally (Euchi et al. 2022;Cinar et al. 2021), mostly due to the increase in life expectancy, as well as the increase in the number of patients with chronic diseases and physical disabilities (Fikar and Hirsch 2017;Tippong et al. 2022). These services not only improve the quality of service provided to the patients but also decrease hospital congestion (Grenouilleau et al. 2019). The pressure on hospitals could be eased by the delivery of a wide range of services from filiation and testing to drug delivery and wound treatment at home. For elderly and disabled patients, vaccination could also be conducted at their homes. While the range of home healthcare services is wide, our problem definition and case study are based on optimizing the routes and schedules of teams that provide filiation services for patients during the pandemic. In our context, filiation services refer to the standard services concerned with assessing the conditions of patients affected by an infectious disease.
In most cases, patients' health conditions differ and it is essential to prioritize them accordingly. To achieve this, a triage level is defined for each patient that shows the degree of urgency and identifies the severity of the patient's condition. These levels are used as weights for the waiting times (latency) of the patients in the objective function of our problem so that urgent patients are more likely to get service sooner. Henceforth, we refer to this routing and scheduling problem as the Home Healthcare Routing and Scheduling Problem with Patient Prioritization (HHRSP-PP).
The HHRSP-PP is closely related to the traveling repairman problem (TRP) and its variants that have been studied in the routing literature. In fact, when the service time of a node (e.g., patient) is added to the traveling time of each incoming arc to that node, the HHRSP-PP can be viewed as a generalization of the multiple traveling repairman problem (mTRP) where the waiting time of each node is multiplied by its priority weight. Recently, the "weighted multiple traveling repairman problem" (WmTRP) has been studied by Muritiba et al. (2021). The authors investigated an application of the WmTRP related to the maintenance of speed cameras, where each node should be visited by only one repairman. In our application, the weight of each node (patient) is associated with the triage level, which is an integer number, typically between 1 and 5.
While the HHRSP-PP and the WmTRP are similar in terms of problem definition, our methodological approach to deal with the problem is different from the approach presented in Muritiba et al. (2021). As opposed to a standard mixed integer programming (MIP) formulation that is only capable of providing exact solutions to very small-sized instances, we develop an alternative level-based integer programming (IP) model, together with valid constraints, by means of which medium-sized instances are solved optimally in reasonable run times. We compare the level-based IP model with the standard MIP model on small-sized instances with 10 and 20 patients (nodes) adopted from the TRP literature (Salehipour et al. 2011).
In order to solve larger instances, we develop a metaheuristic algorithm in which the initial solutions are generated using a problem-specific saving method. A generalized variable neighborhood search (GVNS) procedure is then employed at the improvement step. We test our heuristic algorithm and the level-based IP model on several sets of instances adapted from: (i) the TRP literature (Salehipour et al. 2011), (ii) the vehicle routing problem (VRP) literature (Augerat et al. 1995) and (iii) the WmTRP instances provided in Muritiba et al. (2021). In all the tested instances with up to 60 patients and 15 healthcare service provider teams (HSPTs), optimal solutions are found by both solving the level-based IP model and the proposed heuristic algorithm. However, the metaheuristic algorithm runs within merely 9 s. The performance of our algorithm is further tested on larger instances and compared with well-known algorithms applied to instances of the TRP and the mTRP. We note that our instances and obtained solutions can be used for benchmarking in future studies.
A case study based on the COVID-19 filiation services provided at patient homes in the Kağıthane district of Istanbul is also presented. The data consists of 647 prioritized patients who are clustered into nine regions so that an HSPT is responsible for patients in each cluster. The solution shows that patient prioritization enables patients with more severe conditions to be serviced earlier. A comparison with a centralized optimization approach for a multi-depot system where patients are not clustered prior to optimization reveals how much the centralized system improves both the quality of the solutions and the CPU running times.
The rest of the paper is organized as follows. A review of the related literature is presented in Sect. 2, and the mathematical models are provided in Sect. 3. The metaheuristic algorithm is described in Sect. 4. The results of the experimental study that tests the models and the metaheuristic algorithm are given in Sect. 5. The case study is presented in Sect. 6. Finally, concluding remarks and future work directions are stated in Sect. 7.

Literature review
In this section, we review studies related to home healthcare routing and scheduling, as well as the TRP and its variants.

Home healthcare routing and scheduling
The number of studies addressing home healthcare routing and scheduling has been growing rapidly in recent years. Comprehensive reviews of studies in the context of home healthcare routing and scheduling are provided in Euchi et al. (2022), Grieco et al. (2021), Fikar and Hirsch (2017) and Cisse et al. (2017). Here, we focus on studies that are similar to ours from several perspectives. From the planning horizon viewpoint, the studies on this topic are classified as single-and multi-period. Multi-period models are suitable for problems where a large number of requests are received a priori and must be served on a long-term planning horizon. In this case, the visits must be scheduled over multiple working shifts (Fikar and Hirsch 2017;Bowers et al. 2015). On the other hand, the single-period models are a better fit for problems where all the requests that are accumulated at the beginning of the working shift must be handled within the upcoming shift. In our study, we consider a single-period planning horizon during which all the patients should be visited in the upcoming working shift. Accordingly, the focus of the literature review in this section revolves around single-period studies.
The constraints defined typically in single-period models are time window requirements of the patients, skill requirements, working time regulations and precedence. However, the majority of studies use soft time windows (e.g., Trautsamwieser et al. 2011;Eveborn et al. 2006;Bertels and Fahle 2006) and require all patients to be visited (e.g., Liu et al. (2013); Cappanera et al. (2018)). Given that time windows and latency objectives are in conflict, we do not consider time windows for the patients. Furthermore, we focus on the case with homogeneous teams and consider healthcare service provider teams that are self-contained and can provide all the required services.
The studied problems in the literature can also be categorized as deterministic or stochastic, where stochastic studies address the uncertainty in service times or travel times. Since patients conditions are known a priori in our problem, we assume the service times can be predicted accurately. Thus, our optimization problem has a deterministic nature. As a result, here we focus on deterministic single-period problems. Few studies in this group provide exact solution methods such as branchand-price (e.g., Rasmussen et al. 2012;Manerba and Mansini 2016), and most focus on heuristic approaches. Variable neighborhood search (VNS) (e.g., Mankowska et al. 2014;Trautsamwieser et al. 2011), memetic algorithm (e.g., Decerle et al. 2018), simulated annealing (e.g., Hiermann et al. 2015), genetic algorithm (e.g., Li et al. 2021), particle swarm optimization (e.g., Akjiratikarl et al. 2007), set partitioning heuristics (e.g., Grenouilleau et al. 2019) and quantitative threshold-based approaches (e.g., Nasir and Dang 2020) are among the heuristics developed.
Various objective functions have been considered in the home healthcare routing and scheduling literature, including the minimization of the travel time (e.g., Bredstrom and Ronqvist 2008;Trautsamwieser et al. 2011;Hiermann et al. 2015), travel cost (e.g., Eveborn et al. 2006;Bertels and Fahle 2006;Rasmussen et al. 2012), waiting time of service provider teams (e.g., Trautsamwieser et al. 2011), overtime (e.g., Hiermann et al. 2015;Trautsamwieser et al. 2011), maximization of the preferences (e.g., Bertels and Fahle 2006;Bredstrom and Ronqvist 2008;Trautsamwieser et al. 2011;Rasmussen et al. 2012;Hiermann et al. 2015) and balancing the workload among the HSPTs (e.g., Cappanera et al. 2018;Manerba and Mansini 2016). We depart from these studies by minimizing the weighted total latency of the patients. This type of objective is known to be more challenging than the standard objective of minimizing the total tour length in the routing literature (Angel-Bello et al. 2013). However, in healthcare, waiting can be critical for some patients, and forming routes that minimize the total weighted waiting times carries particular importance in this context. To the best of our knowledge, our article is the first paper that considers the latency objective for a home healthcare routing and scheduling problem.

Traveling repairman problem and its variants
The traveling repairman problem (TRP) aims to keep the customers wait less in a service setting. While the name comes from the repair services, alternative names have also been given by some researchers. In the TRP, a complete graph with a depot node and positive edge distances are provided. The repairman has to perform a route starting from the depot node, visiting all the nodes and ending at the depot node with the objective of minimizing the summation of all visit times. In our problem, the HSPTs must visit all the patients. Moreover, the service time of each patient is incorporated by adding the service time to the traveling time of the incoming arcs to that patient. Since we consider a weighted latency objective function, our problem can be considered as a generalization of the TRP where weights are associated with the customer nodes and multiple repairmen are considered, as opposed to the TRP, which optimizes the route of a single repairman. We next present a brief overview of the previous work on the TRP and its variations.
Numerous articles proposed exact solution procedures including mathematical models for the TRP (e.g., see Sarubbi et al. 2008;Picard and Queyranne 1978;Mendez-Diaz et al. 2008). Mendez-Diaz et al. (2008) formulated a three-index mathematical model that can solve instances having up to 40 nodes by taking advantage of effective valid constraints. Bulhões et al. (2018) developed a branch-and-price algorithm for a set partitioning formulation. Angel-Bello et al. (2013) designed two mathematical models for the TRP based on a multi-level network system. Based on their computational study, Angel-Bello et al. (2013) conclude that level-based models outperform other modeling approaches for the TRP. On the other hand, several studies have developed heuristics and metaheuristics to solve the TRP. Salehipour et al. (2011) developed two metaheuristics where the initial solutions are generated using a greedy randomized adaptive search procedure (GRASP). They have used VNS and variable neighborhood descent (VND) for the improvement phase of their metaheuristics. They called their algorithms the GRASP+VNS and GRASP+VND algorithms. Mladenović et al. (2013) developed a VNS algorithm for the TRP. They tested their algorithm on instances provided in Salehipour et al. (2011) and showed that their proposed algorithm outperforms the GRASP+VNS and GRASP+VND algorithms on these instances.
The multiple TRP (mTRP) is a generalization of the TRP in the sense that it involves finding the routes of multiple repairmen. Luo et al. (2014), Nucamendi et al. (2015, and Nucamendi-Guillén et al. (2016) provided different methods to address the mTRP. For example, Nucamendi-Guillén et al. (2016) proposed a metaheuristic algorithm that contains an iterated greedy (IG) phase for initiation and a local search (LS) phase for improvement. Sze et al. (2017) developed an adaptive variable neighborhood search algorithm to solve the mTRP. Angel-Bello et al. (2019), on the other hand, focused on exact approaches and proposed five mathematical models for the mTRP. The first three formulations emanated from the classical flow-based formulations and the latter two were based on time-dependent models regarding a multi-level network. These models perform much better than the first three. Liu et al. (2018) developed a branch-and-price algorithm to solve the multitrip variation of the mTRP with up to 85 nodes. More recently, larger instances were solved exactly and approximately; Bruni et al. (2022) studied the multi-depot mTRP and proposed two mathematical models as well as a hybrid genetic algorithm to solve instances with up to 240 nodes.
Several authors proposed heuristic algorithms for various variations of the mTRP. Bang (2018) considered a variation of the mTRP with a distance constraint, where each vehicle is not permitted to travel longer than a given limit. They designed a GRASP heuristic to generate the initial solution, and a VND algorithm to improve the initial solution. In Avci and Avci (2019), the authors designed an adaptive large neighborhood search algorithm (ALNS) for solving the mTRP with time-dependent profits. Lalla-Ruiz and Voß (2020) addressed the multi-depot extension of the mTRP and proposed a matheuristic algorithm.
While both the TRP and mTRP are well-studied problems, there are only a few studies that focus on the variations of these problems where weights play a role in the objective function. Dewilde et al. (2013) employed a tabu search algorithm for the prize-collecting TRP, where prizes are associated with the nodes and a subset of the nodes should be selected to be visited to collect the maximum total prize within a given time limit. In this variation, a prize of p i − t i is collected when the repairman arrives at node i ∈ V at time t i , where p i is a fixed prize allocated to node i ∈ V . This variation differs from minimizing the weighted waiting times (latencies) as the problem is selective and the prizes are time-dependent. Lu et al. (2019) developed a memetic algorithm to solve a generalization of the problem introduced in Dewilde et al. (2013) with multiple repairmen. García et al. (2002) provided a linear time exact algorithm and Wu (2000) presented a dynamic programming model for solving a special case of the weighted TRP (WTRP) in which the underlying graph is a path. Akbari and Shiri (2021) addressed an online variation of the weighted TRP considering only one repairman. In a recent study, Muritiba et al. (2021) proposed a branch-and-cut algorithm and developed an iterated local search algorithm to address the weighted multi-repairmen WTRP, motivated by an application involving speed cameras to be maintained.
To the best of our knowledge, none of the TRP variations we have encountered in the literature have been studied in the context of home healthcare. Our paper is the first study that investigates the "latency" objective in this context. That is, we study a TRP with multiple repairmen and a weighted latency objective function, motivated by the home healthcare routing and scheduling problem. Our level-based IP model and metaheuristic algorithm are also new additions to the WmTRP literature.

Mathematical models
We define the HHRSP-PP on a complete graph G = (V 0 , E) with node set V 0 = {0, 1, … , n} and arc set A = {(i, j) ∶ i, j ∈ V 0 , i ≠ j} . The depot (starting location of the HSPTs) is located at node 0 and the set of nodes to be serviced (patients) is represented by V = V 0 ⧵{0} . The distance from node i to node j (or travel time) is given as c ij , (i, j) ∈ A , according to the locations of patients i and j and the shortest path in the road network between them. We incorporate service times in our problem by adding the service time of patient i ( i ∈ {1, 2, … , n}) to the traveling time of each incoming arc to patient i. Since each patient must be visited exactly once, the optimal solution and the optimal objective function value remain unchanged after this transformation. A priority (triage level or weight) is associated with each patient (i.e., each node in V) that is denoted by w i , i ∈ V , where higher values indicate increased priority. In the HHRSP-PP, there are m number of HSPTs that originate from the depot node at time 0, (t 0 = 0) , traverse a number of nodes to service the patients, and return back to the depot. The collection of the routes of the HSPTs should be such that each patient in V is serviced in one of the routes. Clearly, it can be assumed that the number of HSPTs is less than the number of patients that should be serviced, i.e., m < n . Letting t i be the time at which patient i ∈ V is visited, the objective function of the HHRSP-PP is to minimize ∑ i∈V w i t i . Next, we present two formulations for this problem.

Model I
Model I is based on calculating the arrival time to each patient considering the visit time of the previously serviced patients by the same HSPT and the travel time between the patients and also the service time of the considered patient. Since the visit time of a patient must be greater than the time at which the previous patient was visited, the solution cannot contain a closed tour. However, the optimal solution of the problem should include m closed tours starting and ending in the depot node, where m shows the number of HSPTs. In order to remedy this issue, a dummy sink node indexed by n + 1 is defined such that the HSPTs finish their routes at this node. For this purpose, we set the time of going from i ∈ V to the dummy sink node equal to the time of going from i to the depot excluding the service time of node i. We show the union of the node set V together with the sink node n + 1 by V n+1 = V ∪ {n + 1} and the union of all nodes consists of the depot and the sink node by V 0,n+1 . In this model, is the first set of variables that denote whether a HSPT is going from node i to node j or not, and t i ≥ 0 (i ∈ V 0,n+1 ) is the second set of variables that denote the time at which patient i is visited.
The objective function in (1) minimizes the total weighted latency (arrival times) of all the patients nodes. Constraints (2) and (3) ensure that m routes are determined for the m HSPTs. Constraints (4) and (5) ensure that each patient is visited exactly once. Constraints (6) calculate the arrival times and prevent sub-tours. In this constraint, M can get any value greater than or equal to ∑ (i,j)∈A c ij .

Model II
As mentioned, since each patient must be visited exactly once, the HHRSP-PP is equivalent to the WmTRP when the service times of patients are added to the traveling times of incoming arcs to the corresponding patient nodes. Model II is based on a multi-level model developed in Angel-Bello et al. (2019) for the multiple TRP with no weights. We modify the mathematical formulation of Angel-Bello et al.
(2019) to make it applicable to HHRSP-PP. Angel-Bello et al. (2019) showed computationally that a level-based model outperforms flow-based (Gavish and Graves 1978) or multiple TSP-based (Bektas 2006) formulations for latency objectives in terms of computational time. This is because the level-based model is designed to compute the arrival times to the patients (nodes) via the objective function. In Model II, a set of binary variables x r ij are defined where if the arc (i, j) is used to link patient i at level r + 1 with patient j at level r, then it will be equal to one. Levels in this model are defined to facilitate the modeling and show the steps of the movements for each of the HSTPs. An illustration of the level-based model is provided in Fig. 1. In this example, a case with four patients ( n = 4 ) and 2 HSPTs ( m = 2 ) is considered. Since each HSPT visits at least one patient in the optimal solution, the number of levels at which a patient is

Fig. 1
Level-based structure visited is restrained by L = n − m + 1 , which is equal to 4 in our example. A sample route for each of the HSPTs is also presented in Fig. 1. The route of one of the HSPTs include visiting patients 1, 4 and 3 consecutively. For this HSPT, the teams leave the depot in level 4, visits patient v 1 in level 3 and then visits patients v 4 and v 3 in levels 2 and 1, respectively. For the other HSPT, they only meet patient v 2 in level 1.
Given that this level-based formulation was first proposed for the case with no weights on the nodes (patients), to develop Model II we define a transformation of the input graph G = (V 0 , E) to a graph denoted by G = (V 0 , E) . Suppose that the weight (triage level) w i for a node (patient) i in V is an integer number which represents the severity of the condition of patient i. G is a complete graph with node set V 0 that includes the depot and w i duplicated nodes for each node i ∈ V . An example of this transformation is given in Fig. 2. In this example, there are three patients that should be visited. The triage level of the first patient is two (w 1 = 2) , the triage level of the second patient is three ( w 2 = 3 ), and the triage level of the third patient is one (w 3 = 1) . In this transformation of the graph G to graph G , each node i ∈ V with a weight of w i is replaced with w i duplicated nodes denoted by V a i ∶ i ∈ V, a ∈ {1, … , w i } . As a result, G is a complete graph with node set Furthermore, the set of nodes in G excluding the depot is denoted by V.
In graph G , the cost or distance of going from node V a , if both nodes are duplicated nodes of the same patient), then the cost of moving between them equals 0.

Proposition 1
The optimal solution and its objective function value on G is the same as the optimal solution and its objective function value on G.
Proof First, let us assume the optimal visiting times on G are given as t * i for i ∈ V . As a result, ∑ i∈V w i t * i represents the optimal objective function value on G. The dis- . This corresponds to finding a feasible solution on G with the same objective function value obtained on G.
There exists an optimal solution on G such that once one of the w i duplicated nodes of patient i is visited by one HSPT, all the other w i − 1 duplicated nodes of patient i will be visited before that HSPT leaves patient i. This is because the distances between the duplicated nodes for the same patient are set equal to zero and once a HSPT arrives to a duplicated node of a patient, the visiting time of all the other duplicated nodes of that patient is set to that time. This means that an optimal solution on G corresponds to a feasible solution on G with the same objective function value and the proof is completed. ◻ Relying on Proposition 1, we develop a level-based model to solve the HHRSP-PP on G as described in the following. In this model, we only have x r is used by a HSPT to visit node V b j at level r immediately after visiting node V a i at level r + 1 . In this model, the maximum number of levels is L = ∑ i∈V w i − m + 1 . For simplicity of notation, indices u and s are used to show the nodes in set V as well.
The objective function of the level-based model on G is expressed in (9). Constraint (10) guarantees that exactly m HSPTs terminate their routes at level 1. By constraint (11), the number of HSPTs that leave the depot node is set to m. The flow conservation constraints are given by constraints (12) and (13) to ensure the continuity of the paths. Constraint set (14) guarantees that each patient is visited by one and only one HSPT. In the level-based model, this is achieved by ensuring that each patient is active only in one of the levels.

Valid constraints
Since there exist w i duplicated nodes for patient i ( i ∈ {1, 2, … , n} ) in the levelbased formulation, the size of the input graph grows significantly. Moreover, there may be several optimal solutions where the order of the visits of the patients is optimal, but within the duplicated nodes for each patient, the order of visiting the duplicated nodes changes. To remedy this, we propose a number of valid constraints and add them to our level-based model to cut off the feasible space. This way, we eliminate the alternative optimal solutions that correspond to the same order of visiting the patients but allow different combination of visits between the duplicated nodes for the same patient. For that, we set a rule such that for a fixed patient i ∈ V , the HSPT should first visit the duplicated node .., and then V w i i . It should be noted that this rule does not eliminate the optimal solution from the feasible space and only eliminates repetitive solutions that correspond to the same order of visiting the patients. The valid constraints are given in the following.
By (16), the m HSPTs should visit the first indexed duplicated nodes of n patients in V as the first node of their route. Constraints (17) ensure that only the first and last indexed duplicated can be visited first and last of a visit to a patient. Given constraints (18), the HSPTs only visit the dummy nodes allocated to a patient in their indexed order. By constraint (19), there are n − m traversals from the last indexed duplicated node of a patient to the first indexed duplicated node of another patient. By adding constraints (16) to (19) to constraints from (9) to (15), we form the second model named as Model II.

Saving+GVNS heuristic
Model II presented in Sect. 3.2 is shown to solve some moderate-sized instances in our computational tests. This model has two major limitations: firstly, the priority weights (triage levels) of the patients should be integer values and secondly, when the priority weights increase the number of duplicated nodes increases significantly and the problem becomes very difficult to solve. However, here we point out that for our application, the weights of each patient (node) are determined according to triage levels, which are usually integer values between 1 and 5. Nevertheless, in order to overcome these shortcomings, we develop a two-phase algorithm which is capable of handling non-integer weights and is able to solve much larger instances in short time. Similar to Models I and II, we consider the service times in our algorithm by adding the service time of patient i ( i ∈ {1, 2, … , n}) to the traveling time of each incoming arc to the node of patient i in the solutions. At the first phase, a saving procedure is used to generate an initial solution. In the second phase, a general variable neighborhood search (GVNS) algorithm is utilized to improve the initial solution found in the first phase. VNS, which has the advantage that it needs few parameters, has been successfully used in solving diverse problems (see Hansen et al. 2008;Pan et al. 2021;Sadati et al. 2022).
Next, an overview of the proposed algorithm, called the Saving+GVNS algorithm, is presented. The algorithm starts with an initial solution S 0 . A set of neighborhood structures N k (k = 1, … , k max ) are used in the shaking phase and another set of neighborhood structures M l (l = 1, … , l max ) in the local search phase. In the shaking phase, a solution S is generated by applying the first neighborhood N 1 on the initial solution S 0 . For a given solution S , local search is performed by applying the first neighborhood l 1 to obtain a new solution S ′ . If the new solution S ′ has a better objective function value compared to the incumbent solution S * , the incumbent solution S * is updated and local search is continued with the first neighborhood l 1 . Otherwise, l is incremented by 1 and local search is performed using the next neighborhood. This process is continued until all local search neighborhood structures are explored ( l = l max ). In the local search phase, if a new incumbent solution is produced, index k is set to 1. Otherwise, k is incremented by 1 and GVNS restarts from the incumbent solution S * . These steps are repeated until a termination criterion is met. The proposed GVNS is terminated when any of the provided conditions is encountered: (i) a specified number of iterations (10,000 in our experiments), (ii) a number of consecutive non-improving iterations (1000 in our experiments). The pseudo-code of the Saving+GVNS algorithm is given as Algorithm 1. Based on line 9 of Algorithm 1, a solution is feasible if it satisfies constraints (2) to (5) of Model I.

Generating an initial solution using a saving method
In order to construct a feasible initial solution, we apply a modified greedy savings heuristic which is inspired from the Clarke and Wright (CW) savings algorithm (Clarke and Wright 1964). CW is originally proposed for solving the VRP and contains three steps including; (i) initial allocation of one HSPT to each patient, (ii) calculating the saving from merging two routes and sorting the savings in decreasing order and (iii) merging two routes to form a new route starting from the largest saving if it results in a feasible route. Our proposed saving heuristic is different from the original CW in calculations of the saving values and the feasibility of merging two routes. The saving values in our algorithm are calculated as The feasibility of the HHRSP-PP is defined based on using exactly m HSPTs. To achieve this, a pre-determined number of patients are assigned to each HSPT. In this regard, at the merging step of the algorithm, at most ⌈ n m ⌉ patients can be assigned to each HSPT (route) initially, where n is the number of patients and m is the total number of HSPTs. The pseudo-code of the proposed saving procedure for generating initial solutions is presented in Algorithm 2.
In the merging step at line 5 of Algorithm 2, based on the positions of two given patients i and j, four cases can be defined as the following: Case I If patient i is the first and j is the last patient on two separate routes, edges (0, i) and (j, 0) will be removed and a new edge (j, i) will be added to form a single route. Figure 3 gives a visual representation of this merging. Note that in all visual representations of each case (i.e., Figs. 3, 4, 5 and 6), first and second routes give the initial two routes, and the third route shows the merged one.
Case II If patients i and j are both the first patients on their respective routes, initially the order of the first route is reversed and then edges (i, 0) and (0, j) are removed and a new route is formed by adding edge (i, j). Figure 4 gives a visual representation of this merging. Case III If patient i is the last and j is the first patient on their respective routes, edges (i, 0) and (0, j) are eliminated and the merged route is formed by adding (i, j). Figure 5 gives a visual representation of Case III merging.
Case IV If patients i and j are both the last patients on their respective routes, initially the order of the first route is reversed and then edges (0, i) and (j, 0) are removed. Finally, by adding edge (j, i), the merged route is formed. Figure 6 gives a visual representation of this merging.
If the generated number of HSPTs (routes) R at the end of line 11 of Algorithm 2 is greater than m (i.e., m < R ), the extra R-m routes that include the minimum number of patients are discarded from the solution, and the patients in these routes are inserted to the end of the first m routes, which contain the least number of patients (lines 12-14). Note that in line 9 of Algorithm 2, if the new route generated by linking nodes u and v does not satisfy the assignment of at most ⌈ n m ⌉ patients, the route is ignored, and the algorithm continues to check other pairs of nodes in the saving list to form a new route (line 4).

Shaking
In the shaking phase, a solution is generated by applying the set of neighborhoods N k . In our implementation, we used five neighborhoods for the shaking phase of our GVNS. These moves are applied both as intra-route and inter-route moves. In Figs. 7,8,9,10 and 11,parts (a) and (b) show the pre-move and post-move states of the intra-route move, respectively, and parts (c) and (d) show the pre-move and post-move states of the inter-routes move for two HSPTs, respectively. In the case where the routes of two HSPTs are considered (parts (c) and (d) of the figures), the routes are distinguished by solid and dashed lines. 1-0 Move In this move, a patient from his/her current position is removed and inserted in another new position (Fig. 7). 1-1 Exchange This move swaps the positions of two given patients (Fig. 9). 1-2 Move For two given patients, this move swaps the position of the first patient with the second patient and her successor (Fig. 10).
2-2 Exchange For two given patients, this move swaps the positions of the first patient and his successor with the second patient and her successor (Fig. 11).
The aforementioned neighborhoods are implemented both as intra-and interroute moves and are explored in a cyclic sequential order starting with N 1 = 1 -0 move and ending with N 5 = 2 -2 Exchange.

Local search
At each iteration of the GVNS, local search is performed by applying the set of four move operators. In these moves, instead of changing the position of the patients, the steps of the traversed edges are modified.
2-Opt For a given set of two arcs in a single route that define a crisscross, this move substitutes them with two new arcs by reversing the sequence of the nodes visited in between (Fig. 12).
2-Opt* This move is a modification of the 2-Opt move. For a given set of two arcs from two separate routes that create a crisscross, they are replaced with two new edges without reversing the sequence of the patients (Fig. 13). 3-Opt In this move, three arcs in a given route are removed and then the network is reconnected in all feasible ways (Fig. 14).
Or-Opt This move resettles a chain of successive patients by substituting three arcs with three new ones without reversing the sub-routes (Fig. 15).
Here, we point out that 2-Opt, 3-Opt and Or-Opt are applied only as intra-route moves. In these four move operators of local search, we resort to the best improvement approach. These neighborhoods are investigated in a cyclic sequential ordering beginning from M 1 = 2-Opt and terminating with M 4 = Or-Opt. If a new and better solution is produced by applying any of these neighborhood structures, the local search is continued with the first neighborhood M 1 = 2-Opt (i.e., the index l is reinitialized to one). Otherwise, the local search is performed using the next neighborhood. This process is continued until all local search neighborhood structures are explored (lines 7-18 of Algorithm 1).

Computational study
For the purpose of conducting a comprehensive computational experiment, we compare the performance of the proposed Saving+GVNS algorithm with the best known algorithms from the literature for various variants of the WmTRP. All experiments were performed on a computer equipped with Intel Xeon(R) and 3.60 GHz processor. The models were coded in Python and solved using Gurobi Optimizer 9 with an academic license. The Saving+GVNS algorithm was coded in C# Visual Studio 2019.

Performance of the Saving+GVNS algorithm
In this section, we assess the performance of the Saving+GVNS algorithm on simplified versions of WmTRP from the literature. In Sects. 5.1.1 and 5.1.2, our algorithm is compared with TRP and mTRP instances from the literature, respectively.

Performance of the Saving+GVNS on TRP instances
In this section, we analyze the performance of Saving+GVNS algorithm on a TRP variation. This version of the TRP has three main differences with WmTRP: (i) nodes are not weighted and hence the weight of each node is set to one, (ii) only one repairman exists, and (iii) the time in which the repairman returns to the depot is added to the objective function. Moreover, the authors of this study compared their results with the Branch-and-Cut-and-Price (B&C&P) algorithm tested on the TSPLib instances presented in Abeledo et al. (2013). We selected those TSPLib instances that were tested in Mladenović et al. (2013) and utilized two-dimensional Euclidean distances. The results of the experiments are given in Table 1.
In Table 1, the name of the instance is given in the "Instance" column. The objective function values obtained in Abeledo et al. (2013) and Mladenović et al. (2013) were the same in these instances and they are given in the second column. These values are directly extracted from these articles. Those values that are given with a bold font are the optimal solutions of the corresponding instances. The results of the Saving+GVNS algorithm on these instances are given in the next columns. In the presented tables, OFV denotes objective function value and CRT refers to the CPU run time in seconds. Since the VNS of Mladenović et al. (2013) was able to find similar solutions to the B&C&P in all the instances, we only reported the gap between the best found objective function value from the Saving+GVNS algorithm and the results reported in the literature. Since these algorithms were not coded by us and were tested on computers prior to 2013, we have not reported their CPU run time. In total, in nine instances out of 12, optimal solutions were found using the Saving+GVNS algorithm. The maximum gap was in instance rat99 with 1.01% and the average gap over these instances was 0.18%. The average CPU run time over these instances was 132.70 s with a maximum of 246.34 s for the lin105 instance. Considering that the Saving+GVNS algorithm was developed to solve the WmTRP and is now tested on a special case of the WmTRP, obtaining such small gaps is promising and verifies the good performance and robustness of the Saving+GVNS algorithm in terms of both the solution quality and CPU run time.

Performance of the Saving+GVNS on mTRP instances
In this section, we test the performance of the Saving+GVNS algorithm on some mTRP instances that were addressed in Nucamendi-Guillén et al. (2016) and Angel-Bello et al. (2019). The main difference of mTRP and the WmTRP is that in the mTRP, the weight of all the nodes is equal to one. As part of the computational study of these articles, they tested their models and algorithms on 21 P-instances from the VRP literature (Augerat et al. 1995). While Angel-Bello et al. (2019) considered Euclidean distances that were rounded down to the nearest integer, Nucamendi-Guillén et al. (2016) considered Euclidean distances with up to two decimal points. The results of our computational experiments are given in Table 2. The "Instance" column gives the name of the instances, which include the number of nodes and repairmen in that instance. The column denoted by "Angel-Bello et al.
(2019)" gives the extracted values from this article which are obtained for the mTRP using their mathematical model. In the next three columns, the results of testing the Saving+GVNS algorithm on these mTRP instances are given. Our algorithm found the optimal solutions to all of these instances with an average and maximum CRT of 4.84 and just over 11 s, respectively. We note that since Angel-Bello et al. (2019) used a mathematical model, they obtained the optimal solutions to all of these instances. In the next two columns denoted by Nucamendi-Guillén et al. (2016), we present the extracted results of their MIP model and metaheuristic algorithm labeled by IG+LS from their article. Their metaheuristic algorithm finds the optimal solutions to all the tested instances except for P-k8-n23. The results of testing the Saving+GVNS algorithm on these instances are given next. Similar to the case with integer distances, our algorithm found all the optimal solutions to these mTRP instances with an average CRT of under 5 s. While the Saving+GVNS was developed to solve the WmTRP, it is also capable of finding the optimal solutions to all of these mTRP instances verifying its robustness.

WmTRP instances
Next, the Saving+GVNS algorithm and our models are tested on WmTRP instances. This section is divided into two subsections based on the size of the networks and the number of customers in the solved instances.

Small instances
In order to provide a performance comparison of Models I and II, we adopted some small instances from the TRP literature (Salehipour et al. 2011). These small data sets have 10 and 20 nodes. Since these instances are not weighted, based on triage levels that are integer values between 1 and 5, we assigned a uniformly distributed number between 1 and 5 to each of the patient nodes as their weights and then tested both Models I and II and the Saving+GVNS algorithm with 3, 4 and 5 HSPTs on all the instances. The results of testing these instances with 10 and 20 nodes are given in Figs. 16 and 17, respectively. For the instances with 10 nodes, both Models I and II found the optimal solutions in less than 2.5 s in all the instances. The Saving+GVNS algorithm was able to find the same optimal solutions in less than 0.2 s. As a result, Fig. 16 only shows the comparison of the CPU run time between these models and the algorithm. For the case with 20 patients however, Model I was not able to find any of the optimal solutions in a time limit of 1 h. As a result of this, we only gave the CRT comparison between Model II and Saving+GVNS in Fig. 17.
Among the instances with 20 patients, Model II found the optimal solution to all these instances with a maximum CRT of 26.79 s. An observation is that, when the number of HSPTs increases, on the same instance, the problem becomes computationally easier. For example, in instances with 20 nodes, the average CRT over the ten tested instances with 3 HSPTs is 14.51 s, while it reduces to 11.78 and 9.83 s with 4 and 5 HSPTs, respectively. As expected, adding more HSPTs on the same instances reduces the objective function value. While the average objective function value over ten instances with 20 nodes and 3 HSPTs is 3872.9, it reduces to 3408.4 and 3170.5 for 4 and 5 HSPTs, respectively. In both cases with 10 and 20 nodes, the Saving+GVNS algorithm performed found the optimal solution in all the cases. Since the algorithm was able to find the optimal solutions in all the 10 randomized runs for all the instances, the Best OFV and the Avg OFV are the same for all the small instances in this section.

Moderate-sized instances
In this section, we test the performance of the Saving+GVNS and Model II on moderate-sized instances. Some of these instances are directly taken from the literature and some are generated on instances adopted from the VRP literature.

Performance of the Saving+GVNS on WmTRP instances
Recently, Muritiba et al. (2021) introduced the WmTRP and proposed a branch-and-cut (B&C) algorithm to solve this problem. In their study, they addressed the application of the WmTRP for maintenance of speed cameras and introduced six categories of data sets denoted by "brd14051," "d15112," "d18512," "fn14461," "nrw1379" and Fig. 16 Results of the small instances with 10 patients "pr1002." For each of these categories, they generated instances with 30, 40 and 50 nodes and assigned a floating point weight value between 0 and 2 to each of the nodes. For instances with 30 nodes, they used 6 repairmen and for instances with 40 and 50 nodes they used 8 and 10 repairmen, respectively. We note that in our context, these repairmen represent HSPTs and each node represents a patient location. Moreover, for each category and combination of number of nodes and repairmen, they generated 10 instances and reported the average objective function value obtained from their B&C algorithm over them.
The B&C algorithm developed by Muritiba et al. (2021) was able to find the optimal solution in all of their tested instances. Our Saving+GVNS algorithm was also able to find all the same optimal solutions in all the tested instances. Therefore, we only presented the CPU run times in Fig. 18. In this figure, we used the same run times that were presented in Muritiba et al. (2021). For their computational experiments, they coded their B&C algorithm in Java 8 and executed it on an Intel i7-3820 processor with eight cores running at 3.60 GHz and with 16 GB of random access memory, operating under Ubuntu 18.04. As an MILP solver, they used Gurobi 8.0. We can see that the Saving+GVNS algorithm was able to find the optimal solutions in a considerably shorter time. Since Model II only handles integer weights, we have modified the instances described above and generated triage levels for each of the nodes (patients) using integer values between 1 and 5 and solved the new instances using both Model II and our Saving+GVNS algorithm. The results of these experiments are presented in Table 3. Since in Muritiba et al. (2021), the allocated weights to the nodes are floating points between 0 and 2, we modified them to integer weights using w formula where w ′ i shows the weights used by Muritiba et al. (2021) and w i denotes the modified weights that we used in our computational experiments. In order to obtain the optimality gaps of the Saving+GVNS algorithm presented in this table, we calculated them as: Best OFV − Optimal Solution Optimal Solution × 100% . We use the same formula to calculate the optimality gap throughout the article. This term is also referred to as the denominated relative percentage deviation (RPD) in the literature (Vallada  Table 3, all the instances were solved optimally using Model II within a highest CPU time of 20 min. Furthermore, our Saving+GVNS algorithm was able to find the similar optimal solution to all the instances in an average run time of just over 2 s.

VRP instances
In this section, we further analyze the performance of Model II and verify the performance of the Saving+GVNS algorithm. Since Model I was not capable of solving instances with 20 nodes within 1 h time limit, we did not test its performance in this section for the moderate-sized instances. Instead, we adopted 55 instances from the VRP literature (Augerat et al. 1995) with up to 60 nodes and 15 vehicles. We have adopted these instances from three sets of instances denoted by A, B and P in this study. Each instance is named as X-kY-nZ where X denotes whether that instance is from set A, B, or P, Y gives the number of vehicles in that instance and Z gives the number of nodes in the VRP problem. For each of these instances, we used the same number of HSPTs with the number of vehicles in the VRP problem. Furthermore, since these instances are not weighted, we tested them under two different scenarios. In the first scenario, we assigned weights to each node such that the summation of the weights equals 100. This could for example represent triage values under a different scale or a case where importance percentage is allocated to the nodes (patients). In the second scenario, similar to previous parts, we assigned triage levels using uniformly distributed integer values between 1 and 5 for each of the  Tables 4, 5 and 6. For each instance, under each scenario, we report the optimal objective function value and the CPU run time (in seconds) of Model II and denoted them by Optimal OFV and CRT, respectively. For each instance, we also give the results of the Saving+GVNS algorithm including the best objective function value (Best OFV), the CPU run time in seconds and the optimality gap calculated similar to what was explained above. Table 4 gives the results of the experiments on the A-instances. In these instances, the summation of the triage levels for all the patients is assumed to be 100. For the second scenario, however, since the triage level or weight of each patient is a uniformly generated random integer between 1 and 5, the summation of the triage levels over all the patients is variable and hence is given as column W in the tables. As it can be observed in Table 4, Model II solved all the instances of both scenarios optimally, with a maximum CRT of 1302.68 for the A-k9-n55 instance in the second scenario. The maximum CRT of the first scenario belongs to A-k7-n57 with 231.57 s. As it was expected, once the number of nodes and/or the summation of the triage levels over all the patients increases, the CRT of the Model II also increases. This is evident when comparing the same instances over the first and second scenarios. Given that W is fixed at 100 for all the instances in the first scenario, in those second scenario instances that W is more 100, generally, the CPU run time is higher than the first scenario and vice versa. But a comparison between the CRT of the instances shows that in general, the same sized instances with less HSPTs are more difficult to solve. For example, this can be observed comparing instances A-k5-n39 and A-k6-n39. The results given in this table shows the good performance of the Saving+GVNS algorithm. Our metaheuristic algorithm found the optimal solution in all the tested instances and in the largest instance, i.e., A-k9-n60, the optimal solutions of the first and second scenarios were found in 8.23 and 8.10 s, respectively. We note that development of Model II enabled us to find the optimal solutions in these instances which were used to test the performance of the Saving+GVNS algorithm. We recall that the Best OFV values reported in Table 4 are extracted from 10 runs of our Saving+GVNS algorithm over each instance. In some of the instances, the optimal solution was not found in all the 10 runs (i.e., Best OFV is not equal to the Avg OFV). However, the average gap of all the tested A-instances between the Best OFV and the Avg OFV of both scenarios is under 0.01%. Among different instances of the first scenario, instance A-k5-n39 had the highest gap between the Avg OVF and the Best OFV with 0.18%. In the second scenario, instance A-k9-n55 had the largest gap between Avg OFV and Best OFV with 0.08%. Table 5 gives the results of testing Model II and the Saving+GVNS algorithm on B-instances. The columns of this table are the same as those in Table 4. We can see that all the instances were solved optimally using Model II within the 3 h time limit. For instances B-k7-n52 and B-k7-n56 in the second scenario, the CRT of Model II increased to more than 2 h to find the optimal solutions. However, except these two instances and B-k7-n57 of the second scenario, all the remaining instances were solved optimally within 1 h. Among these three instances, B-k7-n52 was solved in less than 62 min. The performance of the Saving+GVNS algorithm is impressive on the B-instances as well. Our metaheuristic algorithm found the optimal solution to all the tested instances with a maximum CRT of just over 8 s for the B-k7-n57 instances. The optimal solutions of the remaining instances were found using our algorithm in less than 8 s. The significant difference between the CRT of Model II and the Saving+GVNS algorithm combined with the fact that our algorithm found the optimal solutions in all of these instances verifies the superior performance of the Saving+GVNS algorithm. In some of the B-instances, not all the 10 runs were able to find the optimal solutions (i.e., Best OFV was not equal to Avg OFV). In the instances of the first and second scenarios, instances B-k7-n50 and B-k6-n43 had the highest gap between Best OFV and Avg OFV with 0.69 and 0.15%, respectively. The average of this gap over all the B-instances remained under 0.05% for the first scenario and under 0.03% for the second scenario. Table 6 gives the results of Model II and the Saving+GVNS algorithm on P-instances. The columns of this table are the same as those of Table 4. The results of the P-instances are similar to those of the A and B instances. Impressively, in all the tested instances, the Saving+GVNS found the optimal solution in less than 8 s. Model II also found the optimal solution to all the tested instances within the given 3 h time limit. Similar to A and B instances, the optimal solution was not found in some runs of a few instances. For the first scenario, the maximum gap between Best OFV and Avg OFV was in instance P-k10-n51 with 0.22% and the average of these gaps was 0.04%. For the second scenario, the average gap between Best OFV and Avg OFV among all the P-instances remained under 0.06% and the maximum gap occurred in instance P-k5-n40 with 0.79%.

Case study
In this section, we provide a case study of our problem related to the filiation operations during the pandemic for the patients residing in the Kağıthane district of the European side of Istanbul. The data consists of 647 patients who are assigned to clusters according to their home addresses, that is, their geographical proximity to each other. The clusters are formed such that each cluster has about the same number of patients and the sum of the distances between the patients in each cluster is minimized by means of a clustering algorithm. Using this clustering method, the 647 patients are clustered into 9 regions. In addition, a priority score based on the triage level which is an integer number between 1 and 5 is assigned to each patient, indicating the urgency of visiting that patient. The available HSPTs (which contain one or more personnel such as a doctor, nurse or other caregivers) are assigned to these regions such that their workloads are balanced and patients that are geographically close to each other are assigned to the same HSPT. Among these nine regions, we selected region 2 with 71 patients and three HSPTs to provide the visualization of the placement of the patients in Figs. 19 and the obtained routes from the Saving+GVNS algorithm for the three HSPTs on the right side of Fig. 20. Here, we note that the traversal time between two coordinates was calculated using the Haversine formula considering an average speed of 30 km/h. For each patient, a service time of 15 min (0.25 h) is considered as their service time. The weighted latency objective function is calculated based on hours each patient have to weight until they are initially met with a HSPT. For the tested instances in this section, we applied our Saving+GVNS algorithm for ten repetitions and reported the best results among them.
The solution obtained by our Saving+GVNS algorithm is depicted on the right side of Fig. 20. The objective function value for visiting these 71 patients with three HSPTs is 512.91. The objective function value is calculated as the total weighted latency of the patients, where the latency of the patients is calculated based on hours. The total traveled distance according to the obtained routes is 16.64 km and the total time including the service times is 18.30 h (i.e., 6.10 h per team on average). As can be observed on the right side of Fig. 20, with a weighted latency  Fig. 20, we gave the illustration of the routes when the weights are ignored, i.e., are set to 1 for all the patients. Although in this solution the total traveled distance reduces to 14.18 km, the corresponding weighted latency objective for these routes is 695.16 which is 35.53% higher than the 512.91 found using a weighted latency objective.
To investigate the impact of the number of HSPTs on the solutions, we solved the region 2 instance using m = 1, 2, 4 and 5 HSPTs as well. The results of these experiments are given in Table 7. In this table, columns "n" and "m" show the number of patients and the number of HSPTs. Columns "OFV", "CRT (s)", "TD (km)", "TT (hr)" and "TT (hr) per team" indicate the objective function value, CPU run time in seconds, total traveled distances, total time (travel and service) and the average total time per team, respectively. According to the obtained results, as expected, once the number of HSPTs increases, the obtained objective function value decreases. We also observe that the rate of this deduction decreases when more HSPTs are considered. For instance, while this deduction is around 34% when instead of one HSPT, two HSPTs are utilized, this deduction is under 18% when the number of HSPTs increases from four to five. Another important observation is based on the last column. While with only one HSPT the tasks will be finished after 18.26 h, when 5 HSPTs are used, the average time of each team decreases to 3.67 h.
In order to investigate the impact of considering patient priorities, we also solved the same instances by considering equal weights for all the patients. In this scenario, we ignored the weights and set the triage levels of all patients equal to 1. This is equivalent to a "un-weighted" problem. We then solved the instance of the second region using different number of HSPTs, i.e., m = 1, 2, 3, 4, and 5. The results are provided in Fig. 21, in which a comparison of the OFV of the un-weighted problem with the original problem (solution provided in Table 7) is given. To be consistent in comparing the results, for the un-weighted cases, we have extracted the routes for each of the HSPTs and calculated the adjusted objective function value using the corresponding triage level for each of the patients. While using the un-weighted latency objective, the routes might be shorter (as depicted in Fig. 21), ignoring the triage level in planning the routes can result in finding inefficient solutions when a weighted latency objective is considered. When the total weighted latency is calculated for these routes, it is seen that the solutions are on average 43.21% worse than the solutions found directly from the weighted objective function. To further analyze the importance of incorporating the triage levels, we provide the average latency (in hours) of patients categorized by triage levels and the number of teams in Table 8. For example, the average latency of patients with w = 1 and m = 1 is 7.68 and 16.54 h in non-weighted and weighted cases, respectively. The results are given in Fig. 22. As it can be observed, for the un-weighted problem, the average latency of the more urgent patients having a triage level of 5 is consistently more than the average latency of the same patients category under a weighted objective. This verifies the importance of considering the triage levels of the patients in this home healthcare delivery problem. However, it can be observed that, in the weighted version, the latency of the patient categories with lower triage levels can be considerably more than the same patients in the un-weighted scenario. For example, when m = 1 , the latency of the least urgent patients with w = 1 , is more than twice in the weighted version when compared to the un-weighted solution.
As mentioned earlier, we selected region 2 among the nine regions to provide its visualization and further analysis. In this part, we present the results on the remaining eight regions using three HSPTs for each of the regions. For each of these regions, a hospital is considered as the depot and hence, for each of these regions, we have a single-depot problem. Moreover, in order to see the effect of a centralized system and to highlight the capability of our proposed GVNS in solving the multi-depot case, we solve the multi-depot centralized case with all patients and the nine depots (hospitals). Similar to the case for each region where three teams were considered, in the multi-depot case we assume three HSPTs are pre-positioned in each of the depots (i.e., in total there are 27 teams). We refer to the first case as "decentralized" and the latter as "centralized." Table 9 provides the results for each of these cases. In this table, the "Region" column indicates the selected region, and the remaining columns are translated in the same way as Table 7. The last row shows the improvements in the centralized cases compared to the decentralized ones. The total objective function value in the decentralized case is 4524.58, whereas in the centralized case it reduces to 4398.96 (a 2.78% improvement). We can also see that solving the centralized instance with 646 patients took 36.48% less time compared to the summation of the times for the decentralized cases over the nine regions. This  analysis shows that considering a centralized system can improve patients' satisfaction and decrease total latency.

Conclusion
We considered a real-life emergency situation such as the COVID-19 outbreak, where standard home healthcare services such as testing or filiation must be provided to the patients at their homes. In this context, we investigated a problem where multiple home healthcare service provider teams (HSPTs) should provide  service to a given set of patients. The patients are prioritized by triage levels based on the severity of their condition or the urgency of servicing them. Each patient should be serviced by exactly one HSPT and the objective is to minimize the summation of the weighted times at which the patients are serviced. We introduced a number of valid constraints for a level-based model driven by a transformation of the input graph and developed an efficient exact model for this problem which is capable of solving moderate-sized instances. We then developed a metaheuristic (Saving+GVNS) algorithm to solve larger instances in shorter time. In our metaheuristic algorithm, the initial solutions are constructed by a problem-specific saving procedure and then improved by a GVNS algorithm.
We performed extensive computational tests. Initially, to compare our Saving+GVNS algorithm with best known algorithms from the literature, we tested its performance on two well-studied special cases of our problem. We then tested our models and the Saving+GVNS algorithm on standard traveling repairman problem and as well as and the vehicle routing problem test instances. To conduct these experiments, we first verified the importance of developing our level-based model to solve even small-sized instances by comparing its performance with a standard MIP model on small-sized instances. We then tested the performance of the Saving+GVNS algorithm by comparing its results with those from a recent article that addresses the weighted K-traveling repairman problem (Muritiba et al. 2021). On these instances, our algorithm was able to find, impressively, all the optimal solutions within merely 2 s.
In order to show the good performance of the Saving+GVNS algorithm, we next compared its results with the level-based model on moderate-sized instances. Our level-based model was able to solve all the moderate-sized instances within a 3 h time limit. Our algorithm was able to find similar optimal solutions to all the small and moderate-sized instances within a maximum CPU run time of 9 s. We also provided a detailed analysis of a real-life case study, based on a data set of a district in the European side of the Istanbul city, which consists of 647 patients. Finally, in order to further analyze the performance of our metaheuristic algorithm, we tested instances with up to 500 nodes and 20 teams and reported the results in Appendix A.
Analyzing the extension of the studied problem where the triage levels of the patients are not known in advance and are only revealed online can be further analyzed. In such cases, not only the triage level will be online, the service time will also be online and it can only be estimated after the conditions of a patient is examined by the HSPTs. Lastly, although we presented our problem and the solution approaches in the context of home health care services, and in particular our case study was related to the filiation services, they can be applied to cases with visits serving other purposes in other application areas, such as repair and maintenance services. The methods can be applied to solve routing and scheduling problems in other check-up services and even distribution problems in city logistics where having a service-oriented objective is important.

Appendix A: Large instances
In this section, we adopted and tested larger instance from the TRP literature (Salehipour et al. 2011) to further analyze the performance of our Saving+GVNS algorithm. These instances were initially generated and tested on the traditional TRP with 1 repairman considering 50, 100, 200 and 500 non-weighted nodes. In order to assign a weight to each of the nodes, which represent triage levels in our study, we generated and assigned a uniformly distributed random number between 1 and 5 to each node. The results of testing the Saving+GVNS algorithm on these instances are given in Tables 10 and 11.
In Table 10, the results of testing the Saving+GVNS algorithm on instances with 50 and 100 nodes are given. In this table, the left hand side columns give the results of the instances with 50 nodes and the columns on the right side give the results of the tested instances with 100 nodes. With both 50 and 100 nodes, 10 instances were introduced in (Salehipour et al. 2011). These instances are distinguished by n − Ri notation, where n shows the number of nodes and i shows the instance number, i ∈ {1, … , 10} . The name of each instance is given in Instance column. The number of HSPTs in each instance is given by m. While with 50 nodes, we tested the Saving+GVNS algorithm with 3, 4 and 5 HSPTs, with 100 nodes, we tested 5, 8 and 10 HSPTs for each instance. We provided the best objective function value and the average objective function value in the next columns. The CPU time of the Saving+GVNS algorithm on average are provided in column CRT. We also presented the average Saving+GVNS gap over the 10 runs of each instance in the Saving+GVNS Gap column. This is calculated as Avg OFV − Best OFV Best OFV × 100%.
As it can be observed, the average CRT of the instances with 50 nodes is 10.91 s and the average Saving+GVNS gap over these instances is 0.22%. The maximum Saving+GVNS gap over instances with 50 nodes is 2.14% for the 50-R3 instance with 3 HSPTs. These low Saving+GVNS gaps verify the robustness of our algorithm. The maximum reported CRT over the instances with 50 nodes is observed in the 50-R8 instance with 15.24 s. We see that the average CRT reduces as the number of HSPTs increases. In an example from instances with 3 HSPTs, the average CRT is 14.05 s, and it lowers to 10.33 and 8.38 s in instances with 4 and 5 HSPTs, respectively.
In the instances with 100 nodes, we set the number of HSPTs equal to 5, 8 and 10. Over these instances, the average CRT is 20.65 s and the average Saving+GVNS gap is 0.49%. As expected, compared to the case with 50 nodes, when the number of nodes increases, both the average CRT and Saving+GVNS gap increases but the average gap remains under 0.49% which again verifies the robustness of our Saving+GVNS algorithm. Over the instances with 100 nodes, the maximum Saving+GVNS gap is 1.68% in the 100-R9 instance with 5 HSPTs. Similar to the instances with 50 nodes, we can observe that the CRT decreases when the number of HSPTs increases. The average CRT over the 10 instances with 5 HSPTs is 29.91 s and it lowers to 17.01 and 15.04 in the instances of 8 and 10 HSPTs, respectively. When the number of HSPTs increases, the achieved best objective function value decreases. However, on average, when the number of HSPTs increases from 5 to 8, the objective function has a decrease of 21.40% and when the number of HSPTs increases from 8 to 10, the average decrease in the objective function is 7.51%. This shows that increasing the number of HSPTs beyond a certain point might not decrease the objective function significantly. Table 11 gives the results of testing the Saving+GVNS algorithm on instances with 200 and 500 nodes. These instances are again adopted from Salehipour et al. (2011) and a uniformly distributed random number between 1 and 5 is assigned to each node as its weight.
The results of the instances with 200 nodes are given on the left hand side of the Table 11. For each of these instances, we set the number of HSPTs 5, 8, 10. The average CRT over these instances is 289.32 s with a maximum of 660.85 s on 200-R5 instance with 5 HSPTs. The average Saving+GVNS gap for the instances with 200 nodes is 1.60%. The average CRT over instances with 5, 8 and 10 HSPTs are 1.84, 1.82 and 1.35%, respectively. This confirms the robustness of the Saving+GVNS algorithm.
For the instances with 500 nodes, we tested each instance with 10, 15 and 20 HSPTs and gave the results on the right side of Table 11. The average CRT over these instances is 1080.66 s. The maximum CRT is up to 1869.74 s for the 500-R10 instance with 10 HSPTs. The average CRT decreases as the number of HSPTs increases. With 10, 15 and 20 HSPTs, the average CRT is 1724.12, 902.37 and 615.48 s, respectively. The average Saving+GVNS gap is 2.17% and this average gap has a decreasing trend as the number of HSPTs increases. While for the instances with 500 nodes and 10 HSPTs, the average Saving+GVNS gap is 3.36%, this average gap decreases to 1.79% for 15 HSPTs and 1.36% with 20 HSPTs.
Data availability In our study, we have used several data sets from the literature. We used TSPLib instances presented in Abeledo et al. (2013) and "P-instances" from Augerat et al. (1995) to show the effectiveness of our algorithm on TRP and mTRP instances. We used the instances from Salehipour et al. (2011) to test the performance of Model I, Model II and our metaheuristic algorithm. We used the instances from Muritiba et al. (2021) to show the effectiveness of our level-based IP model and metaheuristic on a similar problem recently introduced to the literature. We used, modified versions of "A", "B" and "P" instances from Augerat et al. (1995). Finally, we developed our case study data set based on the Kağıthane district in Istanbul, Turkey. This data set includes 647 patients that are clustered into nine groups for meeting their demand. All these data sets are available to be shared with interested researchers upon their request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.